﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using SpiceNet.Circuits;
using SpiceNet.Simulations;
using SpiceNet.Integration;
using SpiceNet.Models.BSIM2p;
using SpiceNet.Diagnostics;
using System.Numerics;
using SpiceNet.Parameters;

namespace SpiceNet.Models
{
    /// <summary>
    /// A BSIM2 mosfet
    /// </summary>
    public class BSIM2 : Instance
    {
        /// <summary>
        /// Get parameters
        /// </summary>
        public override Dictionary<string, IP> ParamTable => _pTable;

        #region Parameters
        private enum _c
        {
            /* device parameters */
            BSIM2_W = 1,
            BSIM2_L = 2,
            BSIM2_AS = 3,
            BSIM2_AD = 4,
            BSIM2_PS = 5,
            BSIM2_PD = 6,
            BSIM2_NRS = 7,
            BSIM2_NRD = 8,
            BSIM2_OFF = 9,
            BSIM2_IC_VBS = 10,
            BSIM2_IC_VDS = 11,
            BSIM2_IC_VGS = 12,
            BSIM2_IC = 13,
            /* device questions */
            BSIM2_DNODE = 241,
            BSIM2_GNODE = 242,
            BSIM2_SNODE = 243,
            BSIM2_BNODE = 244,
            BSIM2_DNODEPRIME = 245,
            BSIM2_SNODEPRIME = 246,
            BSIM2_VBD = 247,
            BSIM2_VBS = 248,
            BSIM2_VGS = 249,
            BSIM2_VDS = 250,
            BSIM2_CD = 251,
            BSIM2_CBS = 252,
            BSIM2_CBD = 253,
            BSIM2_GM = 254,
            BSIM2_GDS = 255,
            BSIM2_GMBS = 256,
            BSIM2_GBD = 257,
            BSIM2_GBS = 258,
            BSIM2_QB = 259,
            BSIM2_CQB = 260,
            BSIM2_QG = 261,
            BSIM2_CQG = 262,
            BSIM2_QD = 263,
            BSIM2_CQD = 264,
            BSIM2_CGG = 265,
            BSIM2_CGD = 266,
            BSIM2_CGS = 267,
            BSIM2_CBG = 268,
            BSIM2_CAPBD = 269,
            BSIM2_CQBD = 270,
            BSIM2_CAPBS = 271,
            BSIM2_CQBS = 272,
            BSIM2_CDG = 273,
            BSIM2_CDD = 274,
            BSIM2_CDS = 275,
            BSIM2_VON = 276,
            BSIM2_QBS = 277,
            BSIM2_QBD = 278,
            BSIM2_SOURCECONDUCT = 279,
            BSIM2_DRAINCONDUCT = 280
        }
        private static Dictionary<string, IP> _pTable = new Dictionary<string, IP>
        {
            { "l", new IP(IF.IOP, (int)_c.BSIM2_L, T.REAL, "Length") },
            { "w", new IP(IF.IOP, (int)_c.BSIM2_W, T.REAL, "Width") },
            { "ad", new IP(IF.IOP, (int)_c.BSIM2_AD, T.REAL, "Drain area") },
            { "as", new IP(IF.IOP, (int)_c.BSIM2_AS, T.REAL, "Source area") },
            { "pd", new IP(IF.IOP, (int)_c.BSIM2_PD, T.REAL, "Drain perimeter") },
            { "ps", new IP(IF.IOP, (int)_c.BSIM2_PS, T.REAL, "Source perimeter") },
            { "nrd", new IP(IF.IOP, (int)_c.BSIM2_NRD, T.REAL, "Number of squares in drain") },
            { "nrs", new IP(IF.IOP, (int)_c.BSIM2_NRS, T.REAL, "Number of squares in source") },
            { "off", new IP(IF.IOP, (int)_c.BSIM2_OFF, T.FLAG, "Device is initially off") },
            { "vds", new IP(IF.IOP, (int)_c.BSIM2_IC_VDS, T.REAL, "Initial D-S voltage") },
            { "vgs", new IP(IF.IOP, (int)_c.BSIM2_IC_VGS, T.REAL, "Initial G-S voltage") },
            { "vbs", new IP(IF.IOP, (int)_c.BSIM2_IC_VBS, T.REAL, "Initial B-S voltage") },
            { "ic", new IP(IF.IP, (int)_c.BSIM2_IC, T.VECTOR, "Vector of DS,GS,BS initial voltages") }
        };

        public Parameter<double> B2l { get; } = new Parameter<double>(5e-6);
        public Parameter<double> B2w { get; } = new Parameter<double>(5e-6);
        public Parameter<double> B2drainArea { get; } = new Parameter<double>();
        public Parameter<double> B2sourceArea { get; } = new Parameter<double>();
        public Parameter<double> B2drainPerimeter { get; } = new Parameter<double>();
        public Parameter<double> B2sourcePerimeter { get; } = new Parameter<double>();
        public Parameter<double> B2drainSquares { get; } = new Parameter<double>(1);
        public Parameter<double> B2sourceSquares { get; } = new Parameter<double>(1);
        public Parameter<double> B2icVDS { get; } = new Parameter<double>();
        public Parameter<double> B2icVGS { get; } = new Parameter<double>();
        public Parameter<double> B2icVBS { get; } = new Parameter<double>();
        public double B2von { get; private set; }
        public double B2vdsat { get; private set; }
        public bool B2off { get; set; }
        #endregion

        #region Constants
        private const int B2vbd = 0;
        private const int B2vbs = 1;
        private const int B2vgs = 2;
        private const int B2vds = 3;
        private const int B2cd = 4;
        private const int B2id = 4;
        private const int B2cbs = 5;
        private const int B2ibs = 5;
        private const int B2cbd = 6;
        private const int B2ibd = 6;
        private const int B2gm = 7;
        private const int B2gds = 8;
        private const int B2gmbs = 9;
        private const int B2gbd = 10;
        private const int B2gbs = 11;
        private const int B2qb = 12;
        private const int B2cqb = 13;
        private const int B2iqb = 13;
        private const int B2qg = 14;
        private const int B2cqg = 15;
        private const int B2iqg = 15;
        private const int B2qd = 16;
        private const int B2cqd = 17;
        private const int B2iqd = 17;
        private const int B2cggb = 18;
        private const int B2cgdb = 19;
        private const int B2cgsb = 20;
        private const int B2cbgb = 21;
        private const int B2cbdb = 22;
        private const int B2cbsb = 23;
        private const int B2capbd = 24;
        private const int B2iqbd = 25;
        private const int B2cqbd = 25;
        private const int B2capbs = 26;
        private const int B2iqbs = 27;
        private const int B2cqbs = 27;
        private const int B2cdgb = 28;
        private const int B2cddb = 29;
        private const int B2cdsb = 30;
        private const int B2vono = 31;
        private const int B2vdsato = 32;
        private const int B2qbs = 33;
        private const int B2qbd = 34;
        #endregion

        /// <summary>
        /// Private variable for this instance
        /// </summary>
        private int B2dNode;  /* number of the gate node of the mosfet */
        private int B2gNode;  /* number of the gate node of the mosfet */
        private int B2sNode;  /* number of the source node of the mosfet */
        private int B2bNode;  /* number of the bulk node of the mosfet */
        private int B2dNodePrime; /* number of the internal drain node of the mosfet */
        private int B2sNodePrime; /* number of the internal source node of the mosfet */
        private int B2state;
        private SizeDependParam param = null;
        private double B2sourceConductance;   /* cond. of source (or 0): set in setup */
        private double B2drainConductance;    /* cond. of drain (or 0): set in setup */
        private int B2mode;

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        public BSIM2(string name) : base(name, 4, typeof(BSIM2Model)) { }

        /// <summary>
        /// Set parameter
        /// </summary>
        protected override void Param(int id, object value, Circuit ckt = null)
        {
            double[] v;
            switch ((_c)id)
            {
                case _c.BSIM2_W:
                    B2w.Par((double)value);
                    break;
                case _c.BSIM2_L:
                    B2l.Par((double)value);
                    break;
                case _c.BSIM2_AS:
                    B2sourceArea.Par((double)value);
                    break;
                case _c.BSIM2_AD:
                    B2drainArea.Par((double)value);
                    break;
                case _c.BSIM2_PS:
                    B2sourcePerimeter.Par((double)value);
                    break;
                case _c.BSIM2_PD:
                    B2drainPerimeter.Par((double)value);
                    break;
                case _c.BSIM2_NRS:
                    B2sourceSquares.Par((double)value);
                    break;
                case _c.BSIM2_NRD:
                    B2drainSquares.Par((double)value);
                    break;
                case _c.BSIM2_OFF:
                    B2off = (bool)value;
                    break;
                case _c.BSIM2_IC_VBS:
                    B2icVBS.Par((double)value);
                    break;
                case _c.BSIM2_IC_VDS:
                    B2icVDS.Par((double)value);
                    break;
                case _c.BSIM2_IC_VGS:
                    B2icVGS.Par((double)value);
                    break;
                case _c.BSIM2_IC:
                    v = (double[])value;
                    switch (v.Length)
                    {
                        case 3:
                            B2icVBS.Par(v[2]); goto case 2;
                        case 2:
                            B2icVGS.Par(v[1]); goto case 1;
                        case 1:
                            B2icVDS.Par(v[0]); break;
                        default:
                            throw new BadParameterException();
                    }
                    break;
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Ask a parameter
        /// </summary>
        protected override object Ask(int id, Circuit ckt = null)
        {
            switch ((_c)id)
            {
                case _c.BSIM2_L:
                    return B2l.Value;
                case _c.BSIM2_W:
                    return B2w.Value;
                case _c.BSIM2_AS:
                    return B2sourceArea.Value;
                case _c.BSIM2_AD:
                    return B2drainArea.Value;
                case _c.BSIM2_PS:
                    return B2sourcePerimeter.Value;
                case _c.BSIM2_PD:
                    return B2drainPerimeter.Value;
                case _c.BSIM2_NRS:
                    return B2sourceSquares.Value;
                case _c.BSIM2_NRD:
                    return B2drainSquares.Value;
                case _c.BSIM2_OFF:
                    return B2off;
                case _c.BSIM2_IC_VBS:
                    return B2icVBS.Value;
                case _c.BSIM2_IC_VDS:
                    return B2icVDS.Value;
                case _c.BSIM2_IC_VGS:
                    return B2icVGS.Value;
                case _c.BSIM2_DNODE:
                    return B2dNode;
                case _c.BSIM2_GNODE:
                    return B2gNode;
                case _c.BSIM2_SNODE:
                    return B2sNode;
                case _c.BSIM2_BNODE:
                    return B2bNode;
                case _c.BSIM2_DNODEPRIME:
                    return B2dNodePrime;
                case _c.BSIM2_SNODEPRIME:
                    return B2sNodePrime;
                case _c.BSIM2_SOURCECONDUCT:
                    return B2sourceConductance;
                case _c.BSIM2_DRAINCONDUCT:
                    return B2drainConductance;
                case _c.BSIM2_VBD:
                    return ckt.State.States[0][B2state + B2vbd];
                case _c.BSIM2_VBS:
                    return ckt.State.States[0][B2state + B2vbs];
                case _c.BSIM2_VGS:
                    return ckt.State.States[0][B2state + B2vgs];
                case _c.BSIM2_VDS:
                    return ckt.State.States[0][B2state + B2vds];
                case _c.BSIM2_CD:
                    return ckt.State.States[0][B2state + B2cd];
                case _c.BSIM2_CBS:
                    return ckt.State.States[0][B2state + B2cbs];
                case _c.BSIM2_CBD:
                    return ckt.State.States[0][B2state + B2cbd];
                case _c.BSIM2_GM:
                    return ckt.State.States[0][B2state + B2gm];
                case _c.BSIM2_GDS:
                    return ckt.State.States[0][B2state + B2gds];
                case _c.BSIM2_GMBS:
                    return ckt.State.States[0][B2state + B2gmbs];
                case _c.BSIM2_GBD:
                    return ckt.State.States[0][B2state + B2gbd];
                case _c.BSIM2_GBS:
                    return ckt.State.States[0][B2state + B2gbs];
                case _c.BSIM2_QB:
                    return ckt.State.States[0][B2state + B2qb];
                case _c.BSIM2_CQB:
                    return ckt.State.States[0][B2state + B2cqb];
                case _c.BSIM2_QG:
                    return ckt.State.States[0][B2state + B2qg];
                case _c.BSIM2_CQG:
                    return ckt.State.States[0][B2state + B2cqg];
                case _c.BSIM2_QD:
                    return ckt.State.States[0][B2state + B2qd];
                case _c.BSIM2_CQD:
                    return ckt.State.States[0][B2state + B2cqd];
                case _c.BSIM2_CGG:
                    return ckt.State.States[0][B2state + B2cggb];
                case _c.BSIM2_CGD:
                    return ckt.State.States[0][B2state + B2cgdb];
                case _c.BSIM2_CGS:
                    return ckt.State.States[0][B2state + B2cgsb];
                case _c.BSIM2_CBG:
                    return ckt.State.States[0][B2state + B2cbgb];
                case _c.BSIM2_CAPBD:
                    return ckt.State.States[0][B2state + B2capbd];
                case _c.BSIM2_CQBD:
                    return ckt.State.States[0][B2state + B2cqbd];
                case _c.BSIM2_CAPBS:
                    return ckt.State.States[0][B2state + B2capbs];
                case _c.BSIM2_CQBS:
                    return ckt.State.States[0][B2state + B2cqbs];
                case _c.BSIM2_CDG:
                    return ckt.State.States[0][B2state + B2cdgb];
                case _c.BSIM2_CDD:
                    return ckt.State.States[0][B2state + B2cddb];
                case _c.BSIM2_CDS:
                    return ckt.State.States[0][B2state + B2cdsb];
                case _c.BSIM2_VON:
                    return ckt.State.States[0][B2state + B2vono];
                case _c.BSIM2_QBS:
                    return ckt.State.States[0][B2state + B2qbs];
                case _c.BSIM2_QBD:
                    return ckt.State.States[0][B2state + B2qbd];
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Setup the BSIM2 transistor
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void Setup(Circuit ckt, Model mod, object p = null)
        {
            BSIM2Model model = mod as BSIM2Model;
            List<NodeType> extra = new List<NodeType>();
            bool addDrain = false;
            if ((model.B2sheetResistance != 0.0) && (B2drainSquares != 0))
            {
                extra.Add(NodeType.Voltage);
                addDrain = true;
            }
            bool addSource = false;
            if ((model.B2sheetResistance != 0.0) && (B2sourceSquares != 0))
            {
                extra.Add(NodeType.Voltage);
                addSource = true;
            }

            var nodes = BindNodes(ckt, extra.ToArray());
            B2dNode = nodes[0].Number;
            B2gNode = nodes[1].Number;
            B2sNode = nodes[2].Number;
            B2bNode = nodes[3].Number;
            int index = 4;
            B2dNodePrime = addDrain ? nodes[index++].Number : B2dNode;
            B2sNodePrime = addSource ? nodes[index++].Number : B2sNode;

            // Allocate states
            B2state = CreateState(ckt, 35);
        }

        /// <summary>
        /// Do temperature-dependent calculations
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void Temperature(Circuit ckt, Model mod, object p = null)
        {
            BSIM2Model model = mod as BSIM2Model;

            double EffectiveLength;
            double EffectiveWidth;
            double CoxWoverL, Inv_L, Inv_W, tmp;

            Tuple<double, double> size_param = new Tuple<double, double>(this.B2w, this.B2l);
            if (model.Params.ContainsKey(size_param))
                this.param = model.Params[size_param];
            else
            {
                param = new SizeDependParam();
                model.Params.Add(size_param, param);

                EffectiveLength = B2l - model.B2deltaL * 1.0e-6;
                EffectiveWidth = B2w - model.B2deltaW * 1.0e-6;

                if (EffectiveLength <= 0)
                    throw new CircuitException(String.Format("B2: mosfet {0}, model {1}: Effective channel length <= 0", Name, model.Name));

                if (EffectiveWidth <= 0)
                    throw new CircuitException(String.Format("B2: mosfet {0}, model {1}: Effective channel width <=0", Name, model.Name));

                Inv_L = 1.0e-6 / EffectiveLength;
                Inv_W = 1.0e-6 / EffectiveWidth;
                param.Width = B2w;
                param.Length = B2l;
                param.B2vfb = model.B2vfb0 + model.B2vfbW * Inv_W
                   + model.B2vfbL * Inv_L;
                param.B2phi = model.B2phi0 + model.B2phiW * Inv_W
                   + model.B2phiL * Inv_L;
                param.B2k1 = model.B2k10 + model.B2k1W * Inv_W
                   + model.B2k1L * Inv_L;
                param.B2k2 = model.B2k20 + model.B2k2W * Inv_W
                   + model.B2k2L * Inv_L;
                param.B2eta0 = model.B2eta00
                   + model.B2eta0W * Inv_W
                   + model.B2eta0L * Inv_L;
                param.B2etaB = model.B2etaB0 + model.B2etaBW
                   * Inv_W + model.B2etaBL * Inv_L;
                param.B2beta0 = model.B2mob00;
                param.B2beta0B = model.B2mob0B0
                   + model.B2mob0BW * Inv_W
                   + model.B2mob0BL * Inv_L;
                param.B2betas0 = model.B2mobs00
                   + model.B2mobs0W * Inv_W
                   + model.B2mobs0L * Inv_L;
                if (param.B2betas0 < 1.01 * param.B2beta0)
                    param.B2betas0 = 1.01 * param.B2beta0;
                param.B2betasB = model.B2mobsB0
                   + model.B2mobsBW * Inv_W
                   + model.B2mobsBL * Inv_L;
                tmp = (param.B2betas0 - param.B2beta0
                       - param.B2beta0B * model.B2vbb);
                if ((-param.B2betasB * model.B2vbb) > tmp)
                    param.B2betasB = -tmp / model.B2vbb;
                param.B2beta20 = model.B2mob200
                  + model.B2mob20W * Inv_W
                  + model.B2mob20L * Inv_L;
                param.B2beta2B = model.B2mob2B0
                  + model.B2mob2BW * Inv_W
                  + model.B2mob2BL * Inv_L;
                param.B2beta2G = model.B2mob2G0
                  + model.B2mob2GW * Inv_W
                  + model.B2mob2GL * Inv_L;
                param.B2beta30 = model.B2mob300
                  + model.B2mob30W * Inv_W
                  + model.B2mob30L * Inv_L;
                param.B2beta3B = model.B2mob3B0
                  + model.B2mob3BW * Inv_W
                  + model.B2mob3BL * Inv_L;
                param.B2beta3G = model.B2mob3G0
                  + model.B2mob3GW * Inv_W
                  + model.B2mob3GL * Inv_L;
                param.B2beta40 = model.B2mob400
                  + model.B2mob40W * Inv_W
                  + model.B2mob40L * Inv_L;
                param.B2beta4B = model.B2mob4B0
                  + model.B2mob4BW * Inv_W
                  + model.B2mob4BL * Inv_L;
                param.B2beta4G = model.B2mob4G0
                  + model.B2mob4GW * Inv_W
                  + model.B2mob4GL * Inv_L;

                CoxWoverL = model.B2Cox * EffectiveWidth / EffectiveLength;

                param.B2beta0 *= CoxWoverL;
                param.B2beta0B *= CoxWoverL;
                param.B2betas0 *= CoxWoverL;
                param.B2betasB *= CoxWoverL;
                param.B2beta30 *= CoxWoverL;
                param.B2beta3B *= CoxWoverL;
                param.B2beta3G *= CoxWoverL;
                param.B2beta40 *= CoxWoverL;
                param.B2beta4B *= CoxWoverL;
                param.B2beta4G *= CoxWoverL;

                param.B2ua0 = model.B2ua00 + model.B2ua0W * Inv_W
                   + model.B2ua0L * Inv_L;
                param.B2uaB = model.B2uaB0 + model.B2uaBW * Inv_W
                   + model.B2uaBL * Inv_L;
                param.B2ub0 = model.B2ub00 + model.B2ub0W * Inv_W
                   + model.B2ub0L * Inv_L;
                param.B2ubB = model.B2ubB0 + model.B2ubBW * Inv_W
                   + model.B2ubBL * Inv_L;
                param.B2u10 = model.B2u100 + model.B2u10W * Inv_W
                   + model.B2u10L * Inv_L;
                param.B2u1B = model.B2u1B0 + model.B2u1BW * Inv_W
                   + model.B2u1BL * Inv_L;
                param.B2u1D = model.B2u1D0 + model.B2u1DW * Inv_W
                   + model.B2u1DL * Inv_L;
                param.B2n0 = model.B2n00 + model.B2n0W * Inv_W
                   + model.B2n0L * Inv_L;
                param.B2nB = model.B2nB0 + model.B2nBW * Inv_W
                   + model.B2nBL * Inv_L;
                param.B2nD = model.B2nD0 + model.B2nDW * Inv_W
                   + model.B2nDL * Inv_L;
                if (param.B2n0 < 0.0)
                    param.B2n0 = 0.0;

                param.B2vof0 = model.B2vof00
                   + model.B2vof0W * Inv_W
                   + model.B2vof0L * Inv_L;
                param.B2vofB = model.B2vofB0
                   + model.B2vofBW * Inv_W
                   + model.B2vofBL * Inv_L;
                param.B2vofD = model.B2vofD0
                   + model.B2vofDW * Inv_W
                   + model.B2vofDL * Inv_L;
                param.B2ai0 = model.B2ai00 + model.B2ai0W * Inv_W
                   + model.B2ai0L * Inv_L;
                param.B2aiB = model.B2aiB0 + model.B2aiBW * Inv_W
                   + model.B2aiBL * Inv_L;
                param.B2bi0 = model.B2bi00 + model.B2bi0W * Inv_W
                   + model.B2bi0L * Inv_L;
                param.B2biB = model.B2biB0 + model.B2biBW * Inv_W
                   + model.B2biBL * Inv_L;
                param.B2vghigh = model.B2vghigh0
                   + model.B2vghighW * Inv_W
                   + model.B2vghighL * Inv_L;
                param.B2vglow = model.B2vglow0
                   + model.B2vglowW * Inv_W
                   + model.B2vglowL * Inv_L;

                param.CoxWL = model.B2Cox * EffectiveLength
                           * EffectiveWidth * 1.0e4;
                param.One_Third_CoxWL = param.CoxWL / 3.0;
                param.Two_Third_CoxWL = 2.0
                           * param.One_Third_CoxWL;
                param.B2GSoverlapCap = model.B2gateSourceOverlapCap
                   * EffectiveWidth;
                param.B2GDoverlapCap = model.B2gateDrainOverlapCap
                   * EffectiveWidth;
                param.B2GBoverlapCap = model.B2gateBulkOverlapCap
                   * EffectiveLength;
                param.SqrtPhi = Math.Sqrt(param.B2phi);
                param.Phis3 = param.SqrtPhi
                               * param.B2phi;
                param.Arg = param.B2betasB
                   - param.B2beta0B - model.B2vdd
                           * (param.B2beta3B - model.B2vdd
                           * param.B2beta4B);
            }


            /* process drain series resistance */
            if ((B2drainConductance = model.B2sheetResistance *
                    B2drainSquares) != 0.0)
            {
                B2drainConductance = 1.0 / B2drainConductance;
            }

            /* process source series resistance */
            if ((B2sourceConductance = model.B2sheetResistance *
                    B2sourceSquares) != 0.0)
            {
                B2sourceConductance = 1.0 / B2sourceConductance;
            }


            param.B2vt0 = param.B2vfb
              + param.B2phi
              + param.B2k1 * param.SqrtPhi
              - param.B2k2 * param.B2phi;
            B2von = param.B2vt0; /* added for initialization*/
        }

        /// <summary>
        /// Load the BSIM2 mosfet
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void Load(Circuit ckt, Model mod, object p = null)
        {
            BSIM2Model model = mod as BSIM2Model;

            double DrainSatCurrent;
            double EffectiveLength;
            double GateBulkOverlapCap;
            double GateDrainOverlapCap;
            double GateSourceOverlapCap;
            double SourceSatCurrent;
            double DrainArea;
            double SourceArea;
            double DrainPerimeter;
            double SourcePerimeter;
            double arg;
            double capbd = 0.0;
            double capbs = 0.0;
            double cbd;
            double cbhat;
            double cbs;
            double cd;
            double cdrain = 0.0;
            double cdhat;
            double cdreq;
            // double ceq;
            double ceqbd;
            double ceqbs;
            double ceqqb;
            double ceqqd;
            double ceqqg;
            double czbd;
            double czbdsw;
            double czbs;
            double czbssw;
            double delvbd;
            double delvbs;
            double delvds;
            double delvgd;
            double delvgs;
            double evbd;
            double evbs;
            double gbd;
            double gbs;
            double gcbdb = 0.0;
            double gcbgb = 0.0;
            double gcbsb = 0.0;
            double gcddb = 0.0;
            double gcdgb = 0.0;
            double gcdsb = 0.0;
            double gcgdb = 0.0;
            double gcggb = 0.0;
            double gcgsb = 0.0;
            double gcsdb = 0.0;
            double gcsgb = 0.0;
            double gcssb = 0.0;
            double gds = 0.0;
            // double geq;
            double gm = 0.0;
            double gmbs = 0.0;
            double sarg;
            double sargsw;
            // double tol;
            double vbd;
            double vbs;
            double vcrit;
            double vds;
            double vdsat;
            double vgb;
            double vgd;
            double vgdo;
            double vgs;
            double von;
            // double xfact;
            double xnrm;
            double xrev;
            bool Check;
            double cgdb = 0.0;
            double cgsb = 0.0;
            double cbdb = 0.0;
            double cdgb = 0.0;
            double cddb = 0.0;
            double cdsb = 0.0;
            double cggb = 0.0;
            double cbgb = 0.0;
            double cbsb = 0.0;
            double csgb = 0.0;
            double cssb = 0.0;
            double csdb = 0.0;
            double PhiB;
            double PhiBSW;
            double MJ;
            double MJSW;
            double argsw;
            double qgate = 0.0;
            double qbulk = 0.0;
            double qdrn = 0.0;
            double qsrc = 0.0;
            double cqgate;
            double cqbulk;
            double cqdrn;
            double vt0;
            double[] args = new double[8];
            CircuitState state = ckt.State;


            EffectiveLength = B2l - model.B2deltaL * 1.0e-6;/* m */
            DrainArea = B2drainArea;
            SourceArea = B2sourceArea;
            DrainPerimeter = B2drainPerimeter;
            SourcePerimeter = B2sourcePerimeter;
            if ((DrainSatCurrent = DrainArea * model.B2jctSatCurDensity)
                    < 1e-15)
            {
                DrainSatCurrent = 1.0e-15;
            }
            if ((SourceSatCurrent = SourceArea * model.B2jctSatCurDensity)
                    < 1.0e-15)
            {
                SourceSatCurrent = 1.0e-15;
            }
            GateSourceOverlapCap = model.B2gateSourceOverlapCap * B2w;
            GateDrainOverlapCap = model.B2gateDrainOverlapCap * B2w;
            GateBulkOverlapCap = model.B2gateBulkOverlapCap * EffectiveLength;
            von = model.B2type * B2von;
            vdsat = model.B2type * B2vdsat;
            vt0 = model.B2type * param.B2vt0;

            Check = true;
            // ByPass = 0;
            if ((ckt.Mode & Circuit.Modes.InitSmSig) != 0)
            {
                vbs = state.States[0][B2state + B2vbs];
                vgs = state.States[0][B2state + B2vgs];
                vds = state.States[0][B2state + B2vds];
            }
            else if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
            {
                vbs = state.States[1][B2state + B2vbs];
                vgs = state.States[1][B2state + B2vgs];
                vds = state.States[1][B2state + B2vds];
            }
            else if ((ckt.Mode & Circuit.Modes.InitJct) != 0 && !B2off)
            {
                vds = model.B2type * B2icVDS;
                vgs = model.B2type * B2icVGS;
                vbs = model.B2type * B2icVBS;
                if ((vds == 0) && (vgs == 0) && (vbs == 0) &&
                        ((ckt.Mode &
                        (Circuit.Modes.Tran | Circuit.Modes.Ac | Circuit.Modes.DcOp | Circuit.Modes.DcTranCurve)) != 0 ||
                        ((ckt.Mode & Circuit.Modes.Uic) == 0)))
                {
                    vbs = -1;
                    vgs = vt0;
                    vds = 0;
                }
            }
            else if ((ckt.Mode & (Circuit.Modes.InitJct | Circuit.Modes.InitFix)) != 0 &&
                  (B2off))
            {
                vbs = vgs = vds = 0;
            }
            else
            {
                vbs = model.B2type * (
                    state.Solution[B2bNode] -
                    state.Solution[B2sNodePrime]);
                vgs = model.B2type * (
                    state.Solution[B2gNode] -
                    state.Solution[B2sNodePrime]);
                vds = model.B2type * (
                    state.Solution[B2dNodePrime] -
                    state.Solution[B2sNodePrime]);
                vbd = vbs - vds;
                vgd = vgs - vds;
                vgdo = state.States[0][B2state + B2vgs] -
                    state.States[0][B2state + B2vds];
                delvbs = vbs - state.States[0][B2state + B2vbs];
                delvbd = vbd - state.States[0][B2state + B2vbd];
                delvgs = vgs - state.States[0][B2state + B2vgs];
                delvds = vds - state.States[0][B2state + B2vds];
                delvgd = vgd - vgdo;

                if (B2mode >= 0)
                {
                    cdhat =
                        state.States[0][B2state + B2cd] -
                        state.States[0][B2state + B2gbd] * delvbd +
                        state.States[0][B2state + B2gmbs] * delvbs +
                        state.States[0][B2state + B2gm] * delvgs +
                        state.States[0][B2state + B2gds] * delvds;
                }
                else
                {
                    cdhat =
                        state.States[0][B2state + B2cd] -
                        (state.States[0][B2state + B2gbd] -
                          state.States[0][B2state + B2gmbs]) * delvbd -
                        state.States[0][B2state + B2gm] * delvgd +
                        state.States[0][B2state + B2gds] * delvds;
                }
                cbhat =
                    state.States[0][B2state + B2cbs] +
                    state.States[0][B2state + B2cbd] +
                    state.States[0][B2state + B2gbd] * delvbd +
                    state.States[0][B2state + B2gbs] * delvbs;

                von = model.B2type * B2von;
                if (state.States[0][B2state + B2vds] >= 0)
                {
                    vgs = fetlim(vgs, state.States[0][B2state + B2vgs]
                            , von);
                    vds = vgs - vgd;
                    vds = limvds(vds, state.States[0][B2state + B2vds]);
                    vgd = vgs - vds;
                }
                else
                {
                    vgd = fetlim(vgd, vgdo, von);
                    vds = vgs - vgd;
                    vds = -limvds(-vds, -(state.States[0][B2state + B2vds]));
                    vgs = vgd + vds;
                }
                if (vds >= 0)
                {
                    vcrit = Circuit.CONSTvt0 * Math.Log(Circuit.CONSTvt0 / (Circuit.CONSTroot2 * SourceSatCurrent));
                    vbs = pnjlim(vbs, state.States[0][B2state + B2vbs],
                            Circuit.CONSTvt0, vcrit, ref Check); /* B2 test */
                    vbd = vbs - vds;
                }
                else
                {
                    vcrit = Circuit.CONSTvt0 * Math.Log(Circuit.CONSTvt0 / (Circuit.CONSTroot2 * DrainSatCurrent));
                    vbd = pnjlim(vbd, state.States[0][B2state + B2vbd],
                            Circuit.CONSTvt0, vcrit, ref Check); /* B2 test*/
                    vbs = vbd + vds;
                }
            }

            /* determine DC current and derivatives */
            vbd = vbs - vds;
            vgd = vgs - vds;
            vgb = vgs - vbs;


            if (vbs <= 0.0)
            {
                gbs = SourceSatCurrent / Circuit.CONSTvt0 + ckt.Config.Gmin;
                cbs = gbs * vbs;
            }
            else
            {
                evbs = Math.Exp(vbs / Circuit.CONSTvt0);
                gbs = SourceSatCurrent * evbs / Circuit.CONSTvt0 + ckt.Config.Gmin;
                cbs = SourceSatCurrent * (evbs - 1) + ckt.Config.Gmin * vbs;
            }
            if (vbd <= 0.0)
            {
                gbd = DrainSatCurrent / Circuit.CONSTvt0 + ckt.Config.Gmin;
                cbd = gbd * vbd;
            }
            else
            {
                evbd = Math.Exp(vbd / Circuit.CONSTvt0);
                gbd = DrainSatCurrent * evbd / Circuit.CONSTvt0 + ckt.Config.Gmin;
                cbd = DrainSatCurrent * (evbd - 1) + ckt.Config.Gmin * vbd;
            }
            /* line 400 */
            if (vds >= 0)
            {
                /* normal mode */
                B2mode = 1;
            }
            else
            {
                /* inverse mode */
                B2mode = -1;
            }
            /* call B2evaluate to calculate drain current and its 
             * derivatives and charge and capacitances related to gate
             * drain, and bulk
             */
            if (vds >= 0)
            {
                B2evaluate(vds, vbs, vgs, ref gm, ref gds, ref gmbs, ref qgate,
                    ref qbulk, ref qdrn, ref cggb, ref cgdb, ref cgsb, ref cbgb, ref cbdb, ref cbsb, ref cdgb,
                    ref cddb, ref cdsb, ref cdrain, ref von, ref vdsat, ckt, model);
            }
            else
            {
                B2evaluate(-vds, vbd, vgd, ref gm, ref gds, ref gmbs, ref qgate,
                    ref qbulk, ref qsrc, ref cggb, ref cgsb, ref cgdb, ref cbgb, ref cbsb, ref cbdb, ref csgb,
                    ref cssb, ref csdb, ref cdrain, ref von, ref vdsat, ckt, model);
            }

            B2von = model.B2type * von;
            B2vdsat = model.B2type * vdsat;



            /*
             *  COMPUTE EQUIVALENT DRAIN CURRENT SOURCE
             */
            cd = B2mode * cdrain - cbd;
            if ((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac | Circuit.Modes.InitSmSig)) != 0 ||
                    ((ckt.Mode & Circuit.Modes.TranOp) != 0 &&
                    (ckt.Mode & Circuit.Modes.Uic) != 0))
            {
                /*
                 *  charge storage elements
                 *
                 *   bulk-drain and bulk-source depletion capacitances
                 *  czbd : zero bias drain junction capacitance
                 *  czbs : zero bias source junction capacitance
                 * czbdsw:zero bias drain junction sidewall capacitance
                 * czbssw:zero bias source junction sidewall capacitance
                 */

                czbd = model.B2unitAreaJctCap * DrainArea;
                czbs = model.B2unitAreaJctCap * SourceArea;
                czbdsw = model.B2unitLengthSidewallJctCap * DrainPerimeter;
                czbssw = model.B2unitLengthSidewallJctCap * SourcePerimeter;
                PhiB = model.B2bulkJctPotential;
                PhiBSW = model.B2sidewallJctPotential;
                MJ = model.B2bulkJctBotGradingCoeff;
                MJSW = model.B2bulkJctSideGradingCoeff;

                /* Source Bulk Junction */
                if (vbs < 0)
                {
                    arg = 1 - vbs / PhiB;
                    argsw = 1 - vbs / PhiBSW;
                    sarg = Math.Exp(-MJ * Math.Log(arg));
                    sargsw = Math.Exp(-MJSW * Math.Log(argsw));
                    state.States[0][B2state + B2qbs] =
                        PhiB * czbs * (1 - arg * sarg) / (1 - MJ) + PhiBSW *
                    czbssw * (1 - argsw * sargsw) / (1 - MJSW);
                    capbs = czbs * sarg + czbssw * sargsw;
                }
                else
                {
                    state.States[0][B2state + B2qbs] =
                        vbs * (czbs + czbssw) + vbs * vbs * (czbs * MJ * 0.5 / PhiB
                        + czbssw * MJSW * 0.5 / PhiBSW);
                    capbs = czbs + czbssw + vbs * (czbs * MJ / PhiB +
                        czbssw * MJSW / PhiBSW);
                }

                /* Drain Bulk Junction */
                if (vbd < 0)
                {
                    arg = 1 - vbd / PhiB;
                    argsw = 1 - vbd / PhiBSW;
                    sarg = Math.Exp(-MJ * Math.Log(arg));
                    sargsw = Math.Exp(-MJSW * Math.Log(argsw));
                    state.States[0][B2state + B2qbd] =
                        PhiB * czbd * (1 - arg * sarg) / (1 - MJ) + PhiBSW *
                    czbdsw * (1 - argsw * sargsw) / (1 - MJSW);
                    capbd = czbd * sarg + czbdsw * sargsw;
                }
                else
                {
                    state.States[0][B2state + B2qbd] =
                        vbd * (czbd + czbdsw) + vbd * vbd * (czbd * MJ * 0.5 / PhiB
                        + czbdsw * MJSW * 0.5 / PhiBSW);
                    capbd = czbd + czbdsw + vbd * (czbd * MJ / PhiB +
                        czbdsw * MJSW / PhiBSW);
                }

            }




            /*
             *  check convergence
             */
            if ((!B2off) || ((ckt.Mode & Circuit.Modes.InitFix) == 0))
            {
                if (Check)
                    ckt.IsCon = false;
            }
            state.States[0][B2state + B2vbs] = vbs;
            state.States[0][B2state + B2vbd] = vbd;
            state.States[0][B2state + B2vgs] = vgs;
            state.States[0][B2state + B2vds] = vds;
            state.States[0][B2state + B2cd] = cd;
            state.States[0][B2state + B2cbs] = cbs;
            state.States[0][B2state + B2cbd] = cbd;
            state.States[0][B2state + B2gm] = gm;
            state.States[0][B2state + B2gds] = gds;
            state.States[0][B2state + B2gmbs] = gmbs;
            state.States[0][B2state + B2gbd] = gbd;
            state.States[0][B2state + B2gbs] = gbs;

            state.States[0][B2state + B2cggb] = cggb;
            state.States[0][B2state + B2cgdb] = cgdb;
            state.States[0][B2state + B2cgsb] = cgsb;

            state.States[0][B2state + B2cbgb] = cbgb;
            state.States[0][B2state + B2cbdb] = cbdb;
            state.States[0][B2state + B2cbsb] = cbsb;

            state.States[0][B2state + B2cdgb] = cdgb;
            state.States[0][B2state + B2cddb] = cddb;
            state.States[0][B2state + B2cdsb] = cdsb;

            state.States[0][B2state + B2capbs] = capbs;
            state.States[0][B2state + B2capbd] = capbd;

            /* bulk and channel charge plus overlaps */

            if (((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac)) == 0) &&
                    (((ckt.Mode & Circuit.Modes.TranOp) == 0) ||
                    ((ckt.Mode & Circuit.Modes.Uic) == 0)) && ((ckt.Mode
                    & Circuit.Modes.InitSmSig) == 0)) goto line850;

        // line755:
            if (B2mode > 0)
            {

                args[0] = GateDrainOverlapCap;
                args[1] = GateSourceOverlapCap;
                args[2] = GateBulkOverlapCap;
                args[3] = capbd;
                args[4] = capbs;
                args[5] = cggb;
                args[6] = cgdb;
                args[7] = cgsb;

                B2mosCap(ckt, vgd, vgs, vgb,
            args,
                    /*
                    GateDrainOverlapCap,
                    GateSourceOverlapCap,GateBulkOverlapCap,
                    capbd,capbs,cggb,cgdb,cgsb,
                    */
            cbgb, cbdb, cbsb, cdgb, cddb, cdsb
            , ref gcggb, ref gcgdb, ref gcgsb, ref gcbgb, ref gcbdb, ref gcbsb, ref gcdgb
            , ref gcddb, ref gcdsb, ref gcsgb, ref gcsdb, ref gcssb, ref qgate, ref qbulk
            , ref qdrn, ref qsrc);
            }
            else
            {

                args[0] = GateSourceOverlapCap;
                args[1] = GateDrainOverlapCap;
                args[2] = GateBulkOverlapCap;
                args[3] = capbs;
                args[4] = capbd;
                args[5] = cggb;
                args[6] = cgsb;
                args[7] = cgdb;

                B2mosCap(ckt, vgs, vgd, vgb, args,
                    /*
                    GateSourceOverlapCap,
                            GateDrainOverlapCap,GateBulkOverlapCap,
                    capbs,capbd,cggb,cgsb,cgdb,
                    */
            cbgb, cbsb, cbdb, csgb, cssb, csdb
                    , ref gcggb, ref gcgsb, ref gcgdb, ref gcbgb, ref gcbsb, ref gcbdb, ref gcsgb
                    , ref gcssb, ref gcsdb, ref gcdgb, ref gcdsb, ref gcddb, ref qgate, ref qbulk
                    , ref qsrc, ref qdrn);
            }

            // if (ByPass) goto line860;
            state.States[0][B2state + B2qg] = qgate;
            state.States[0][B2state + B2qd] = qdrn -
                    state.States[0][B2state + B2qbd];
            state.States[0][B2state + B2qb] = qbulk +
                    state.States[0][B2state + B2qbd] +
                    state.States[0][B2state + B2qbs];

            /* store small signal parameters */
            if (((ckt.Mode & (Circuit.Modes.Ac | Circuit.Modes.Tran)) == 0) &&
                    (ckt.Mode & Circuit.Modes.TranOp) != 0 && (ckt.Mode &
                    Circuit.Modes.Uic) != 0) goto line850;
            if ((ckt.Mode & Circuit.Modes.InitSmSig) != 0)
            {
                state.States[0][B2state + B2cggb] = cggb;
                state.States[0][B2state + B2cgdb] = cgdb;
                state.States[0][B2state + B2cgsb] = cgsb;
                state.States[0][B2state + B2cbgb] = cbgb;
                state.States[0][B2state + B2cbdb] = cbdb;
                state.States[0][B2state + B2cbsb] = cbsb;
                state.States[0][B2state + B2cdgb] = cdgb;
                state.States[0][B2state + B2cddb] = cddb;
                state.States[0][B2state + B2cdsb] = cdsb;
                state.States[0][B2state + B2capbd] = capbd;
                state.States[0][B2state + B2capbs] = capbs;

                goto line1000;
            }

            if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
            {
                state.States[1][B2state + B2qb] =
                    state.States[0][B2state + B2qb];
                state.States[1][B2state + B2qg] =
                    state.States[0][B2state + B2qg];
                state.States[1][B2state + B2qd] =
                    state.States[0][B2state + B2qd];
            }

            ckt.Integration.Integrate(state, B2state + B2qb, 0.0);
            ckt.Integration.Integrate(state, B2state + B2qg, 0.0);
            ckt.Integration.Integrate(state, B2state + B2qd, 0.0);

            goto line860;

        line850:
            /* initialize to zero charge conductance and current */
            ceqqg = ceqqb = ceqqd = 0.0;
            gcdgb = gcddb = gcdsb = 0.0;
            gcsgb = gcsdb = gcssb = 0.0;
            gcggb = gcgdb = gcgsb = 0.0;
            gcbgb = gcbdb = gcbsb = 0.0;
            goto line900;

        line860:
            /* evaluate equivalent charge current */
            cqgate = state.States[0][B2state + B2iqg];
            cqbulk = state.States[0][B2state + B2iqb];
            cqdrn = state.States[0][B2state + B2iqd];
            ceqqg = cqgate - gcggb * vgb + gcgdb * vbd + gcgsb * vbs;
            ceqqb = cqbulk - gcbgb * vgb + gcbdb * vbd + gcbsb * vbs;
            ceqqd = cqdrn - gcdgb * vgb + gcddb * vbd + gcdsb * vbs;

            if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
            {
                state.States[1][B2state + B2iqb] =
                    state.States[0][B2state + B2iqb];
                state.States[1][B2state + B2iqg] =
                    state.States[0][B2state + B2iqg];
                state.States[1][B2state + B2iqd] =
                    state.States[0][B2state + B2iqd];
            }

            /*
             *  load current vector
             */
        line900:

            ceqbs = model.B2type * (cbs - (gbs - ckt.Config.Gmin) * vbs);
            ceqbd = model.B2type * (cbd - (gbd - ckt.Config.Gmin) * vbd);

            ceqqg = model.B2type * ceqqg;
            ceqqb = model.B2type * ceqqb;
            ceqqd = model.B2type * ceqqd;
            if (B2mode >= 0)
            {
                xnrm = 1;
                xrev = 0;
                cdreq = model.B2type * (cdrain - gds * vds - gm * vgs - gmbs * vbs);
            }
            else
            {
                xnrm = 0;
                xrev = 1;
                cdreq = -(model.B2type) * (cdrain + gds * vds - gm * vgd - gmbs * vbd);
            }

            state.Rhs[B2gNode] -= ceqqg;
            state.Rhs[B2bNode] -= (ceqbs + ceqbd + ceqqb);
            state.Rhs[B2dNodePrime] +=
                    (ceqbd - cdreq - ceqqd);
            state.Rhs[B2sNodePrime] +=
                    (cdreq + ceqbs + ceqqg + ceqqb + ceqqd);

            /*
             *  load y matrix
             */

            state.Matrix[B2dNode, B2dNode] += (B2drainConductance);
            state.Matrix[B2gNode, B2gNode] += (gcggb);
            state.Matrix[B2sNode, B2sNode] += (B2sourceConductance);
            state.Matrix[B2bNode, B2bNode] += (gbd + gbs - gcbgb - gcbdb - gcbsb);
            state.Matrix[B2dNodePrime, B2dNodePrime] +=
                (B2drainConductance + gds + gbd + xrev * (gm + gmbs) + gcddb);
            state.Matrix[B2sNodePrime, B2sNodePrime] +=
                (B2sourceConductance + gds + gbs + xnrm * (gm + gmbs) + gcssb);
            state.Matrix[B2dNode, B2dNodePrime] += (-B2drainConductance);
            state.Matrix[B2gNode, B2bNode] += (-gcggb - gcgdb - gcgsb);
            state.Matrix[B2gNode, B2dNodePrime] += (gcgdb);
            state.Matrix[B2gNode, B2sNodePrime] += (gcgsb);
            state.Matrix[B2sNode, B2sNodePrime] += (-B2sourceConductance);
            state.Matrix[B2bNode, B2gNode] += (gcbgb);
            state.Matrix[B2bNode, B2dNodePrime] += (-gbd + gcbdb);
            state.Matrix[B2bNode, B2sNodePrime] += (-gbs + gcbsb);
            state.Matrix[B2dNodePrime, B2dNode] += (-B2drainConductance);
            state.Matrix[B2dNodePrime, B2gNode] += ((xnrm - xrev) * gm + gcdgb);
            state.Matrix[B2dNodePrime, B2bNode] += (-gbd + (xnrm - xrev) * gmbs - gcdgb - gcddb - gcdsb);
            state.Matrix[B2dNodePrime, B2sNodePrime] += (-gds - xnrm * (gm + gmbs) + gcdsb);
            state.Matrix[B2sNodePrime, B2gNode] += (-(xnrm - xrev) * gm + gcsgb);
            state.Matrix[B2sNodePrime, B2sNode] += (-B2sourceConductance);
            state.Matrix[B2sNodePrime, B2bNode] += (-gbs - (xnrm - xrev) * gmbs - gcsgb - gcsdb - gcssb);
            state.Matrix[B2sNodePrime, B2dNodePrime] += (-gds - xrev * (gm + gmbs) + gcsdb);


        line1000: ;

        }

        /// <summary>
        /// B2 Evaluate
        /// </summary>
        private void B2evaluate(double Vds, double Vbs, double Vgs, ref double gm, ref double gds, ref double gmb, ref double qg, ref double qb, ref double qd, ref double cgg, ref double cgd, ref double cgs,
            ref double cbg, ref double cbd, ref double cbs, ref double cdg, ref double cdd, ref double cds, ref double Ids, ref double von, ref double vdsat, Circuit ckt, BSIM2Model model)
        {
            double Vth, Vdsat = 0.0;
            double Phisb, T1s, Eta, Gg, Aa, Inv_Aa, U1, U1s, Vc, Kk, SqrtKk;
            double dPhisb_dVb, dT1s_dVb, dVth_dVb, dVth_dVd, dAa_dVb, dVc_dVd;
            double dVc_dVg, dVc_dVb, dKk_dVc, dVdsat_dVd = 0.0, dVdsat_dVg = 0.0, dVdsat_dVb = 0.0;
            double dUvert_dVg, dUvert_dVd, dUvert_dVb, Inv_Kk; // dBeta0_dVd, dBeta0_dVb, 
            double dUtot_dVd, dUtot_dVb, dUtot_dVg, Ai, Bi, Vghigh, Vglow, Vgeff, Vof;
            double Vbseff, Vgst, Vgdt, Qbulk, Utot;
            double T0, T1, T2, T3, T4, T5, Arg1, Arg2, Exp0 = 0.0;
            double tmp, tmp1, tmp2, tmp3, Uvert, Beta1, Beta2, Beta0, dGg_dVb, Exp1 = 0.0;
            double T6, T7, T8, T9, n = 0.0, ExpArg, ExpArg1; // Usat
            double Beta, dQbulk_dVb, dVgdt_dVg, dVgdt_dVd;
            double dVbseff_dVb, Ua, Ub, dVgdt_dVb, dQbulk_dVd;
            double Con1, Con3, Con4, SqrVghigh, SqrVglow, CubVghigh, CubVglow;
            double delta, Coeffa, Coeffb, Coeffc, Coeffd, Inv_Uvert, Inv_Utot;
            double Inv_Vdsat, tanh, Sqrsech, dBeta1_dVb, dU1_dVd, dU1_dVg, dU1_dVb;
            double Betaeff, FR, dFR_dVd, dFR_dVg, dFR_dVb, Betas, Beta3, Beta4;
            double dBeta_dVd, dBeta_dVg, dBeta_dVb, dVgeff_dVg, dVgeff_dVd, dVgeff_dVb;
            double dCon3_dVd, dCon3_dVb, dCon4_dVd, dCon4_dVb, dCoeffa_dVd, dCoeffa_dVb;
            double dCoeffb_dVd, dCoeffb_dVb, dCoeffc_dVd, dCoeffc_dVb;
            double dCoeffd_dVd, dCoeffd_dVb;
            bool ChargeComputationNeeded;
            int valuetypeflag;			/* added  3/19/90 JSD   */

            if ((ckt.Mode & (Circuit.Modes.Ac | Circuit.Modes.Tran)) != 0 ||
                ((ckt.Mode & Circuit.Modes.TranOp) != 0 && (ckt.Mode & Circuit.Modes.Uic) != 0) ||
                (ckt.Mode & Circuit.Modes.InitSmSig) != 0)
            {
                ChargeComputationNeeded = true;
            }
            else
            {
                ChargeComputationNeeded = false;
            }

            if (Vbs < model.B2vbb2) Vbs = model.B2vbb2;
            if (Vgs > model.B2vgg2) Vgs = model.B2vgg2;
            if (Vds > model.B2vdd2) Vds = model.B2vdd2;

            /* Threshold Voltage. */
            if (Vbs <= 0.0)
            {
                Phisb = param.B2phi - Vbs;
                dPhisb_dVb = -1.0;
                T1s = Math.Sqrt(Phisb);
                dT1s_dVb = -0.5 / T1s;
            }
            else
            {
                tmp = param.B2phi / (param.B2phi + Vbs);
                Phisb = param.B2phi * tmp;
                dPhisb_dVb = -tmp * tmp;
                T1s = param.Phis3 / (param.B2phi + 0.5 * Vbs);
                dT1s_dVb = -0.5 * T1s * T1s / param.Phis3;
            }

            Eta = param.B2eta0 + param.B2etaB * Vbs;
            Ua = param.B2ua0 + param.B2uaB * Vbs;
            Ub = param.B2ub0 + param.B2ubB * Vbs;
            U1s = param.B2u10 + param.B2u1B * Vbs;

            Vth = param.B2vfb + param.B2phi + param.B2k1
                * T1s - param.B2k2 * Phisb - Eta * Vds;
            dVth_dVd = -Eta;
            dVth_dVb = param.B2k1 * dT1s_dVb + param.B2k2
                 - param.B2etaB * Vds;

            Vgst = Vgs - Vth;

            tmp = 1.0 / (1.744 + 0.8364 * Phisb);
            Gg = 1.0 - tmp;
            dGg_dVb = 0.8364 * tmp * tmp * dPhisb_dVb;
            T0 = Gg / T1s;
            tmp1 = 0.5 * T0 * param.B2k1;
            Aa = 1.0 + tmp1;
            dAa_dVb = (Aa - 1.0) * (dGg_dVb / Gg - dT1s_dVb / T1s);
            Inv_Aa = 1.0 / Aa;

            Vghigh = param.B2vghigh;
            Vglow = param.B2vglow;

            if ((Vgst >= Vghigh) || (param.B2n0 == 0.0))
            {
                Vgeff = Vgst;
                dVgeff_dVg = 1.0;
                dVgeff_dVd = -dVth_dVd;
                dVgeff_dVb = -dVth_dVb;
            }
            else
            {
                Vof = param.B2vof0 + param.B2vofB * Vbs
                + param.B2vofD * Vds;
                n = param.B2n0 + param.B2nB / T1s
                  + param.B2nD * Vds;
                tmp = 0.5 / (n * model.B2Vtm);

                ExpArg1 = -Vds / model.B2Vtm;
                ExpArg1 = Math.Max(ExpArg1, -30.0);
                Exp1 = Math.Exp(ExpArg1);
                tmp1 = 1.0 - Exp1;
                tmp1 = Math.Max(tmp1, 1.0e-18);
                tmp2 = 2.0 * Aa * tmp1;

                if (Vgst <= Vglow)
                {
                    ExpArg = Vgst * tmp;
                    ExpArg = Math.Max(ExpArg, -30.0);
                    Exp0 = Math.Exp(0.5 * Vof + ExpArg);
                    Vgeff = Math.Sqrt(tmp2) * model.B2Vtm * Exp0;
                    T0 = n * model.B2Vtm;
                    dVgeff_dVg = Vgeff * tmp;
                    dVgeff_dVd = dVgeff_dVg * (n / tmp1 * Exp1 - dVth_dVd - Vgst
                           * param.B2nD / n + T0 * param.B2vofD);
                    dVgeff_dVb = dVgeff_dVg * (param.B2vofB * T0
                           - dVth_dVb + param.B2nB * Vgst
                           / (n * T1s * T1s) * dT1s_dVb + T0 * Inv_Aa * dAa_dVb);
                }
                else
                {
                    ExpArg = Vglow * tmp;
                    ExpArg = Math.Max(ExpArg, -30.0);
                    Exp0 = Math.Exp(0.5 * Vof + ExpArg);
                    Vgeff = Math.Sqrt(2.0 * Aa * (1.0 - Exp1)) * model.B2Vtm * Exp0;
                    Con1 = Vghigh;
                    Con3 = Vgeff;
                    Con4 = Con3 * tmp;
                    SqrVghigh = Vghigh * Vghigh;
                    SqrVglow = Vglow * Vglow;
                    CubVghigh = Vghigh * SqrVghigh;
                    CubVglow = Vglow * SqrVglow;
                    T0 = 2.0 * Vghigh;
                    T1 = 2.0 * Vglow;
                    T2 = 3.0 * SqrVghigh;
                    T3 = 3.0 * SqrVglow;
                    T4 = Vghigh - Vglow;
                    T5 = SqrVghigh - SqrVglow;
                    T6 = CubVghigh - CubVglow;
                    T7 = Con1 - Con3;
                    delta = (T1 - T0) * T6 + (T2 - T3) * T5 + (T0 * T3 - T1 * T2) * T4;
                    delta = 1.0 / delta;
                    Coeffb = (T1 - Con4 * T0) * T6 + (Con4 * T2 - T3) * T5
                   + (T0 * T3 - T1 * T2) * T7;
                    Coeffc = (Con4 - 1.0) * T6 + (T2 - T3) * T7 + (T3 - Con4 * T2) * T4;
                    Coeffd = (T1 - T0) * T7 + (1.0 - Con4) * T5 + (Con4 * T0 - T1) * T4;
                    Coeffa = SqrVghigh * (Coeffc + Coeffd * T0);
                    Vgeff = (Coeffa + Vgst * (Coeffb + Vgst * (Coeffc + Vgst * Coeffd)))
                  * delta;
                    dVgeff_dVg = (Coeffb + Vgst * (2.0 * Coeffc + 3.0 * Vgst * Coeffd))
                           * delta;
                    T7 = Con3 * tmp;
                    T8 = dT1s_dVb * param.B2nB / (T1s * T1s * n);
                    T9 = n * model.B2Vtm;
                    dCon3_dVd = T7 * (n * Exp1 / tmp1 - Vglow * param.B2nD
                          / n + T9 * param.B2vofD);
                    dCon3_dVb = T7 * (T9 * Inv_Aa * dAa_dVb + Vglow * T8
                          + T9 * param.B2vofB);
                    dCon4_dVd = tmp * dCon3_dVd - T7 * param.B2nD / n;
                    dCon4_dVb = tmp * dCon3_dVb + T7 * T8;

                    dCoeffb_dVd = dCon4_dVd * (T2 * T5 - T0 * T6) + dCon3_dVd
                    * (T1 * T2 - T0 * T3);
                    dCoeffc_dVd = dCon4_dVd * (T6 - T2 * T4) + dCon3_dVd * (T3 - T2);
                    dCoeffd_dVd = dCon4_dVd * (T0 * T4 - T5) + dCon3_dVd * (T0 - T1);
                    dCoeffa_dVd = SqrVghigh * (dCoeffc_dVd + dCoeffd_dVd * T0);

                    dVgeff_dVd = -dVgeff_dVg * dVth_dVd + (dCoeffa_dVd + Vgst
                           * (dCoeffb_dVd + Vgst * (dCoeffc_dVd + Vgst
                           * dCoeffd_dVd))) * delta;

                    dCoeffb_dVb = dCon4_dVb * (T2 * T5 - T0 * T6) + dCon3_dVb
                    * (T1 * T2 - T0 * T3);
                    dCoeffc_dVb = dCon4_dVb * (T6 - T2 * T4) + dCon3_dVb * (T3 - T2);
                    dCoeffd_dVb = dCon4_dVb * (T0 * T4 - T5) + dCon3_dVb * (T0 - T1);
                    dCoeffa_dVb = SqrVghigh * (dCoeffc_dVb + dCoeffd_dVb * T0);

                    dVgeff_dVb = -dVgeff_dVg * dVth_dVb + (dCoeffa_dVb + Vgst
                           * (dCoeffb_dVb + Vgst * (dCoeffc_dVb + Vgst
                           * dCoeffd_dVb))) * delta;
                }
            }

            if (Vgeff > 0.0)
            {
                Uvert = 1.0 + Vgeff * (Ua + Vgeff * Ub);
                Uvert = Math.Max(Uvert, 0.2);
                Inv_Uvert = 1.0 / Uvert;
                T8 = Ua + 2.0 * Ub * Vgeff;
                dUvert_dVg = T8 * dVgeff_dVg;
                dUvert_dVd = T8 * dVgeff_dVd;
                dUvert_dVb = T8 * dVgeff_dVb + Vgeff * (param.B2uaB
                       + Vgeff * param.B2ubB);

                T8 = U1s * Inv_Aa * Inv_Uvert;
                Vc = T8 * Vgeff;
                T9 = Vc * Inv_Uvert;
                dVc_dVg = T8 * dVgeff_dVg - T9 * dUvert_dVg;
                dVc_dVd = T8 * dVgeff_dVd - T9 * dUvert_dVd;
                dVc_dVb = T8 * dVgeff_dVb + param.B2u1B * Vgeff * Inv_Aa
                    * Inv_Uvert - Vc * Inv_Aa * dAa_dVb - T9 * dUvert_dVb;


                tmp2 = Math.Sqrt(1.0 + 2.0 * Vc);
                Kk = 0.5 * (1.0 + Vc + tmp2);
                Inv_Kk = 1.0 / Kk;
                dKk_dVc = 0.5 + 0.5 / tmp2;
                SqrtKk = Math.Sqrt(Kk);

                T8 = Inv_Aa / SqrtKk;
                Vdsat = Vgeff * T8;
                Vdsat = Math.Max(Vdsat, 1.0e-18);
                Inv_Vdsat = 1.0 / Vdsat;
                T9 = 0.5 * Vdsat * Inv_Kk * dKk_dVc;
                dVdsat_dVd = T8 * dVgeff_dVd - T9 * dVc_dVd;
                dVdsat_dVg = T8 * dVgeff_dVg - T9 * dVc_dVg;
                dVdsat_dVb = T8 * dVgeff_dVb - T9 * dVc_dVb - Vdsat * Inv_Aa * dAa_dVb;

                Beta0 = param.B2beta0 + param.B2beta0B * Vbs;
                Betas = param.B2betas0 + param.B2betasB * Vbs;
                Beta2 = param.B2beta20 + param.B2beta2B * Vbs
                  + param.B2beta2G * Vgs;
                Beta3 = param.B2beta30 + param.B2beta3B * Vbs
                  + param.B2beta3G * Vgs;
                Beta4 = param.B2beta40 + param.B2beta4B * Vbs
                  + param.B2beta4G * Vgs;
                Beta1 = Betas - (Beta0 + model.B2vdd * (Beta3 - model.B2vdd
                  * Beta4));

                T0 = Vds * Beta2 * Inv_Vdsat;
                T0 = Math.Min(T0, 30.0);
                T1 = Math.Exp(T0);
                T2 = T1 * T1;
                T3 = T2 + 1.0;
                tanh = (T2 - 1.0) / T3;
                Sqrsech = 4.0 * T2 / (T3 * T3);

                Beta = Beta0 + Beta1 * tanh + Vds * (Beta3 - Beta4 * Vds);
                T4 = Beta1 * Sqrsech * Inv_Vdsat;
                T5 = model.B2vdd * tanh;
                dBeta_dVd = Beta3 - 2.0 * Beta4 * Vds + T4 * (Beta2 - T0 * dVdsat_dVd);
                dBeta_dVg = T4 * (param.B2beta2G * Vds - T0 * dVdsat_dVg)
                      + param.B2beta3G * (Vds - T5)
                  - param.B2beta4G * (Vds * Vds - model.B2vdd * T5);
                dBeta1_dVb = param.Arg;
                dBeta_dVb = param.B2beta0B + dBeta1_dVb * tanh + Vds
                      * (param.B2beta3B - Vds * param.B2beta4B)
                      + T4 * (param.B2beta2B * Vds - T0 * dVdsat_dVb);


                if (Vgst > Vglow)
                {
                    if (Vds <= Vdsat) /* triode region */
                    {
                        T3 = Vds * Inv_Vdsat;
                        T4 = T3 - 1.0;
                        T2 = 1.0 - param.B2u1D * T4 * T4;
                        U1 = U1s * T2;
                        Utot = Uvert + U1 * Vds;
                        Utot = Math.Max(Utot, 0.5);
                        Inv_Utot = 1.0 / Utot;
                        T5 = 2.0 * U1s * param.B2u1D * Inv_Vdsat * T4;
                        dU1_dVd = T5 * (T3 * dVdsat_dVd - 1.0);
                        dU1_dVg = T5 * T3 * dVdsat_dVg;
                        dU1_dVb = T5 * T3 * dVdsat_dVb + param.B2u1B * T2;
                        dUtot_dVd = dUvert_dVd + U1 + Vds * dU1_dVd;
                        dUtot_dVg = dUvert_dVg + Vds * dU1_dVg;
                        dUtot_dVb = dUvert_dVb + Vds * dU1_dVb;

                        tmp1 = (Vgeff - 0.5 * Aa * Vds);
                        tmp3 = tmp1 * Vds;
                        Betaeff = Beta * Inv_Utot;
                        Ids = Betaeff * tmp3;
                        T6 = Ids / Betaeff * Inv_Utot;
                        gds = T6 * (dBeta_dVd - Betaeff * dUtot_dVd) + Betaeff * (tmp1
                            + (dVgeff_dVd - 0.5 * Aa) * Vds);
                        gm = T6 * (dBeta_dVg - Betaeff * dUtot_dVg) + Betaeff * Vds * dVgeff_dVg;
                        gmb = T6 * (dBeta_dVb - Betaeff * dUtot_dVb) + Betaeff * Vds
                            * (dVgeff_dVb - 0.5 * Vds * dAa_dVb);
                    }
                    else  /* Saturation */
                    {
                        tmp1 = Vgeff * Inv_Aa * Inv_Kk;
                        tmp3 = 0.5 * Vgeff * tmp1;
                        Betaeff = Beta * Inv_Uvert;
                        Ids = Betaeff * tmp3;
                        T0 = Ids / Betaeff * Inv_Uvert;
                        T1 = Betaeff * Vgeff * Inv_Aa * Inv_Kk;
                        T2 = Ids * Inv_Kk * dKk_dVc;

                        if (param.B2ai0 != 0.0)
                        {
                            Ai = param.B2ai0 + param.B2aiB * Vbs;
                            Bi = param.B2bi0 + param.B2biB * Vbs;
                            T5 = Bi / (Vds - Vdsat);
                            T5 = Math.Min(T5, 30.0);
                            T6 = Math.Exp(-T5);
                            FR = 1.0 + Ai * T6;
                            T7 = T5 / (Vds - Vdsat);
                            T8 = (1.0 - FR) * T7;
                            dFR_dVd = T8 * (dVdsat_dVd - 1.0);
                            dFR_dVg = T8 * dVdsat_dVg;
                            dFR_dVb = T8 * dVdsat_dVb + T6 * (param.B2aiB - Ai
                                * param.B2biB / (Vds - Vdsat));

                            gds = (T0 * (dBeta_dVd - Betaeff * dUvert_dVd) + T1
                             * dVgeff_dVd - T2 * dVc_dVd) * FR + Ids * dFR_dVd;
                            gm = (T0 * (dBeta_dVg - Betaeff * dUvert_dVg)
                            + T1 * dVgeff_dVg - T2 * dVc_dVg) * FR + Ids * dFR_dVg;
                            gmb = (T0 * (dBeta_dVb - Betaeff * dUvert_dVb) + T1
                             * dVgeff_dVb - T2 * dVc_dVb - Ids * Inv_Aa * dAa_dVb)
                             * FR + Ids * dFR_dVb;
                            Ids *= FR;
                        }
                        else
                        {
                            gds = T0 * (dBeta_dVd - Betaeff * dUvert_dVd) + T1
                             * dVgeff_dVd - T2 * dVc_dVd;
                            gm = T0 * (dBeta_dVg - Betaeff * dUvert_dVg) + T1 * dVgeff_dVg
                            - T2 * dVc_dVg;
                            gmb = T0 * (dBeta_dVb - Betaeff * dUvert_dVb) + T1
                             * dVgeff_dVb - T2 * dVc_dVb - Ids * Inv_Aa * dAa_dVb;
                        }
                    } /* end of Saturation */
                }
                else
                {
                    T0 = Exp0 * Exp0;
                    T1 = Exp1;
                    Ids = Beta * model.B2Vtm * model.B2Vtm * T0 * (1.0 - T1);
                    T2 = Ids / Beta;
                    T4 = n * model.B2Vtm;
                    T3 = Ids / T4;
                    if ((Vds > Vdsat) && param.B2ai0 != 0.0)
                    {
                        Ai = param.B2ai0 + param.B2aiB * Vbs;
                        Bi = param.B2bi0 + param.B2biB * Vbs;
                        T5 = Bi / (Vds - Vdsat);
                        T5 = Math.Min(T5, 30.0);
                        T6 = Math.Exp(-T5);
                        FR = 1.0 + Ai * T6;
                        T7 = T5 / (Vds - Vdsat);
                        T8 = (1.0 - FR) * T7;
                        dFR_dVd = T8 * (dVdsat_dVd - 1.0);
                        dFR_dVg = T8 * dVdsat_dVg;
                        dFR_dVb = T8 * dVdsat_dVb + T6 * (param.B2aiB - Ai
                            * param.B2biB / (Vds - Vdsat));
                    }
                    else
                    {
                        FR = 1.0;
                        dFR_dVd = 0.0;
                        dFR_dVg = 0.0;
                        dFR_dVb = 0.0;
                    }

                    gds = (T2 * dBeta_dVd + T3 * (param.B2vofD * T4 - dVth_dVd
                     - param.B2nD * Vgst / n) + Beta * model.B2Vtm
                 * T0 * T1) * FR + Ids * dFR_dVd;
                    gm = (T2 * dBeta_dVg + T3) * FR + Ids * dFR_dVg;
                    gmb = (T2 * dBeta_dVb + T3 * (param.B2vofB * T4 - dVth_dVb
                     + param.B2nB * Vgst / (n * T1s * T1s) * dT1s_dVb)) * FR
                     + Ids * dFR_dVb;
                    Ids *= FR;
                }
            }
            else
            {
                Ids = 0.0;
                gm = 0.0;
                gds = 0.0;
                gmb = 0.0;
            }


        // ChargeComputation:

            /* Some Limiting of DC Parameters */
            gds = Math.Max(gds, 1.0e-20);


            if ((model.B2channelChargePartitionFlag > 1)
             || ((!ChargeComputationNeeded) &&
             (model.B2channelChargePartitionFlag > -5)))
            {
                qg = 0.0;
                qd = 0.0;
                qb = 0.0;
                cgg = 0.0;
                cgs = 0.0;
                cgd = 0.0;
                cdg = 0.0;
                cds = 0.0;
                cdd = 0.0;
                cbg = 0.0;
                cbs = 0.0;
                cbd = 0.0;
                goto finished;
            }
            else
            {
                if (Vbs < 0.0)
                {
                    Vbseff = Vbs;
                    dVbseff_dVb = 1.0;
                }
                else
                {
                    Vbseff = param.B2phi - Phisb;
                    dVbseff_dVb = -dPhisb_dVb;
                }
                Arg1 = Vgs - Vbseff - param.B2vfb;
                Arg2 = Arg1 - Vgst;
                Qbulk = param.One_Third_CoxWL * Arg2;
                dQbulk_dVb = param.One_Third_CoxWL * (dVth_dVb - dVbseff_dVb);
                dQbulk_dVd = param.One_Third_CoxWL * dVth_dVd;
                if (Arg1 <= 0.0)
                {
                    qg = param.CoxWL * Arg1;
                    qb = -(qg);
                    qd = 0.0;

                    cgg = param.CoxWL;
                    cgd = 0.0;
                    cgs = -cgg * (1.0 - dVbseff_dVb);

                    cdg = 0.0;
                    cdd = 0.0;
                    cds = 0.0;

                    cbg = -param.CoxWL;
                    cbd = 0.0;
                    cbs = -cgs;
                }
                else if (Vgst <= 0.0)
                {
                    T2 = Arg1 / Arg2;
                    T3 = T2 * T2 * (param.CoxWL - param.Two_Third_CoxWL
                   * T2);

                    qg = param.CoxWL * Arg1 * (1.0 - T2 * (1.0 - T2 / 3.0));
                    qb = -(qg);
                    qd = 0.0;

                    cgg = param.CoxWL * (1.0 - T2 * (2.0 - T2));
                    tmp = T3 * dVth_dVb - (cgg + T3) * dVbseff_dVb;
                    cgd = T3 * dVth_dVd;
                    cgs = -(cgg + cgd + tmp);

                    cdg = 0.0;
                    cdd = 0.0;
                    cds = 0.0;

                    cbg = -cgg;
                    cbd = -cgd;
                    cbs = -cgs;
                }
                else
                {
                    if (Vgst < param.B2vghigh)
                    {
                        Uvert = 1.0 + Vgst * (Ua + Vgst * Ub);
                        Uvert = Math.Max(Uvert, 0.2);
                        Inv_Uvert = 1.0 / Uvert;
                        dUvert_dVg = Ua + 2.0 * Ub * Vgst;
                        dUvert_dVd = -dUvert_dVg * dVth_dVd;
                        dUvert_dVb = -dUvert_dVg * dVth_dVb + Vgst
                     * (param.B2uaB + Vgst * param.B2ubB);

                        T8 = U1s * Inv_Aa * Inv_Uvert;
                        Vc = T8 * Vgst;
                        T9 = Vc * Inv_Uvert;
                        dVc_dVg = T8 - T9 * dUvert_dVg;
                        dVc_dVd = -T8 * dVth_dVd - T9 * dUvert_dVd;
                        dVc_dVb = -T8 * dVth_dVb + param.B2u1B * Vgst * Inv_Aa
                            * Inv_Uvert - Vc * Inv_Aa * dAa_dVb - T9 * dUvert_dVb;

                        tmp2 = Math.Sqrt(1.0 + 2.0 * Vc);
                        Kk = 0.5 * (1.0 + Vc + tmp2);
                        Inv_Kk = 1.0 / Kk;
                        dKk_dVc = 0.5 + 0.5 / tmp2;
                        SqrtKk = Math.Sqrt(Kk);

                        T8 = Inv_Aa / SqrtKk;
                        Vdsat = Vgst * T8;
                        T9 = 0.5 * Vdsat * Inv_Kk * dKk_dVc;
                        dVdsat_dVd = -T8 * dVth_dVd - T9 * dVc_dVd;
                        dVdsat_dVg = T8 - T9 * dVc_dVg;
                        dVdsat_dVb = -T8 * dVth_dVb - T9 * dVc_dVb
                       - Vdsat * Inv_Aa * dAa_dVb;
                    }
                    if (Vds >= Vdsat)
                    {       /* saturation region */
                        cgg = param.Two_Third_CoxWL;
                        cgd = -cgg * dVth_dVd + dQbulk_dVd;
                        tmp = -cgg * dVth_dVb + dQbulk_dVb;
                        cgs = -(cgg + cgd + tmp);

                        cbg = 0.0;
                        cbd = -dQbulk_dVd;
                        cbs = dQbulk_dVd + dQbulk_dVb;

                        cdg = -0.4 * cgg;
                        tmp = -cdg * dVth_dVb;
                        cdd = -cdg * dVth_dVd;
                        cds = -(cdg + cdd + tmp);

                        qb = -Qbulk;
                        qg = param.Two_Third_CoxWL * Vgst + Qbulk;
                        qd = cdg * Vgst;
                    }
                    else
                    {       /* linear region  */
                        T7 = Vds / Vdsat;
                        T8 = Vgst / Vdsat;
                        T6 = T7 * T8;
                        T9 = 1.0 - T7;
                        Vgdt = Vgst * T9;
                        T0 = Vgst / (Vgst + Vgdt);
                        T1 = Vgdt / (Vgst + Vgdt);
                        T5 = T0 * T1;
                        T2 = 1.0 - T1 + T5;
                        T3 = 1.0 - T0 + T5;

                        dVgdt_dVg = T9 + T6 * dVdsat_dVg;
                        dVgdt_dVd = T6 * dVdsat_dVd - T8 - T9 * dVth_dVd;
                        dVgdt_dVb = T6 * dVdsat_dVb - T9 * dVth_dVb;

                        qg = param.Two_Third_CoxWL * (Vgst + Vgdt
                           - Vgdt * T0) + Qbulk;
                        qb = -Qbulk;
                        qd = -param.One_Third_CoxWL * (0.2 * Vgdt
                       + 0.8 * Vgst + Vgdt * T1
                       + 0.2 * T5 * (Vgdt - Vgst));

                        cgg = param.Two_Third_CoxWL * (T2 + T3 * dVgdt_dVg);
                        tmp = dQbulk_dVb + param.Two_Third_CoxWL * (T3 * dVgdt_dVb
                            - T2 * dVth_dVb);
                        cgd = param.Two_Third_CoxWL * (T3 * dVgdt_dVd
                         - T2 * dVth_dVd) + dQbulk_dVd;
                        cgs = -(cgg + cgd + tmp);

                        T2 = 0.8 - 0.4 * T1 * (2.0 * T1 + T0 + T0 * (T1 - T0));
                        T3 = 0.2 + T1 + T0 * (1.0 - 0.4 * T0 * (T1 + 3.0 * T0));
                        cdg = -param.One_Third_CoxWL * (T2 + T3 * dVgdt_dVg);
                        tmp = param.One_Third_CoxWL * (T2 * dVth_dVb
                        - T3 * dVgdt_dVb);
                        cdd = param.One_Third_CoxWL * (T2 * dVth_dVd
                         - T3 * dVgdt_dVd);
                        cds = -(cdg + tmp + cdd);

                        cbg = 0.0;
                        cbd = -dQbulk_dVd;
                        cbs = dQbulk_dVd + dQbulk_dVb;
                    }
                }
            }

        finished:       /* returning Values to Calling Routine */
            valuetypeflag = (int)model.B2channelChargePartitionFlag;
            switch (valuetypeflag)
            {
                case 0: Ids = Math.Max(Ids, 1e-50);
                    break;
                case -1: Ids = Math.Max(Ids, 1e-50);
                    break;
                case -2: Ids = gm;
                    break;
                case -3: Ids = gds;
                    break;
                case -4: Ids = 1.0 / gds;
                    break;
                case -5: Ids = gmb;
                    break;
                case -6: Ids = qg / 1.0e-12;
                    break;
                case -7: Ids = qb / 1.0e-12;
                    break;
                case -8: Ids = qd / 1.0e-12;
                    break;
                case -9: Ids = -(qb + qg + qd) / 1.0e-12;
                    break;
                case -10: Ids = cgg / 1.0e-12;
                    break;
                case -11: Ids = cgd / 1.0e-12;
                    break;
                case -12: Ids = cgs / 1.0e-12;
                    break;
                case -13: Ids = -(cgg + cgd + cgs) / 1.0e-12;
                    break;
                case -14: Ids = cbg / 1.0e-12;
                    break;
                case -15: Ids = cbd / 1.0e-12;
                    break;
                case -16: Ids = cbs / 1.0e-12;
                    break;
                case -17: Ids = -(cbg + cbd + cbs) / 1.0e-12;
                    break;
                case -18: Ids = cdg / 1.0e-12;
                    break;
                case -19: Ids = cdd / 1.0e-12;
                    break;
                case -20: Ids = cds / 1.0e-12;
                    break;
                case -21: Ids = -(cdg + cdd + cds) / 1.0e-12;
                    break;
                case -22: Ids = -(cgg + cdg + cbg) / 1.0e-12;
                    break;
                case -23: Ids = -(cgd + cdd + cbd) / 1.0e-12;
                    break;
                case -24: Ids = -(cgs + cds + cbs) / 1.0e-12;
                    break;
                default: Ids = Math.Max(Ids, 1.0e-50);
                    break;
            }
            von = Vth;
            vdsat = Vdsat;
        }

        /// <summary>
        /// B2 mosCap
        /// </summary>
        private void B2mosCap(Circuit ckt, double vgd, double vgs, double vgb,
            double[] args, double cbgb, double cbdb, double cbsb, double cdgb, double cddb, double cdsb,
            ref double gcggbPointer, ref double gcgdbPointer, ref double gcgsbPointer, ref double gcbgbPointer, ref double gcbdbPointer,
            ref double gcbsbPointer, ref double gcdgbPointer, ref double gcddbPointer, ref double gcdsbPointer,
            ref double gcsgbPointer, ref double gcsdbPointer, ref double gcssbPointer, ref double qGatePointer, ref double qBulkPointer,
            ref double qDrainPointer, ref double qSourcePointer)
        {
            double qgd;
            double qgs;
            double qgb;
            double ag0;

            ag0 = ckt.Integration.Delta;
            /* compute equivalent conductance */
            gcdgbPointer = (cdgb - args[0]) * ag0;
            gcddbPointer = (cddb + args[3] + args[0]) * ag0;
            gcdsbPointer = cdsb * ag0;
            gcsgbPointer = -(args[5] + cbgb + cdgb + args[1]) * ag0;
            gcsdbPointer = -(args[6] + cbdb + cddb) * ag0;
            gcssbPointer = (args[4] + args[1] -
                (args[7] + cbsb + cdsb)) * ag0;
            gcggbPointer = (args[5] + args[0] +
                args[1] + args[2]) * ag0;
            gcgdbPointer = (args[6] - args[0]) * ag0;
            gcgsbPointer = (args[7] - args[1]) * ag0;
            gcbgbPointer = (cbgb - args[2]) * ag0;
            gcbdbPointer = (cbdb - args[3]) * ag0;
            gcbsbPointer = (cbsb - args[4]) * ag0;

            /* compute total terminal charge */
            qgd = args[0] * vgd;
            qgs = args[1] * vgs;
            qgb = args[2] * vgb;
            qGatePointer = qGatePointer + qgd + qgs + qgb;
            qBulkPointer = qBulkPointer - qgb;
            qDrainPointer = qDrainPointer - qgd;
            qSourcePointer = -(qGatePointer + qBulkPointer + qDrainPointer);
        }


        /// <summary>
        /// Load small-signal components
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void AcLoad(Circuit ckt, Model mod, object p = null)
        {
            BSIM2Model model = mod as BSIM2Model;

            CircuitStateComplex state = ckt.StateComplex;
            CircuitState cstate = ckt.State;

            int xnrm;
            int xrev;
            double gdpr;
            double gspr;
            double gm;
            double gds;
            double gmbs;
            double gbd;
            double gbs;
            double capbd;
            double capbs;
            double xcggb;
            double xcgdb;
            double xcgsb;
            double xcbgb;
            double xcbdb;
            double xcbsb;
            double xcddb;
            double xcssb;
            double xcdgb;
            double xcsgb;
            double xcdsb;
            double xcsdb;
            double cggb;
            double cgdb;
            double cgsb;
            double cbgb;
            double cbdb;
            double cbsb;
            double cddb;
            double cdgb;
            double cdsb;

            if (B2mode >= 0)
            {
                xnrm = 1;
                xrev = 0;
            }
            else
            {
                xnrm = 0;
                xrev = 1;
            }
            gdpr = B2drainConductance;
            gspr = B2sourceConductance;
            gm = cstate.States[0][B2state + B2gm];
            gds = cstate.States[0][B2state + B2gds];
            gmbs = cstate.States[0][B2state + B2gmbs];
            gbd = cstate.States[0][B2state + B2gbd];
            gbs = cstate.States[0][B2state + B2gbs];
            capbd = cstate.States[0][B2state + B2capbd];
            capbs = cstate.States[0][B2state + B2capbs];
            /*
             *    charge oriented model parameters
             */

            cggb = cstate.States[0][B2state + B2cggb];
            cgsb = cstate.States[0][B2state + B2cgsb];
            cgdb = cstate.States[0][B2state + B2cgdb];

            cbgb = cstate.States[0][B2state + B2cbgb];
            cbsb = cstate.States[0][B2state + B2cbsb];
            cbdb = cstate.States[0][B2state + B2cbdb];

            cdgb = cstate.States[0][B2state + B2cdgb];
            cdsb = cstate.States[0][B2state + B2cdsb];
            cddb = cstate.States[0][B2state + B2cddb];

            xcdgb = (cdgb - param.B2GDoverlapCap);
            xcddb = (cddb + capbd + param.B2GDoverlapCap);
            xcdsb = cdsb;
            xcsgb = -(cggb + cbgb + cdgb + param.B2GSoverlapCap);
            xcsdb = -(cgdb + cbdb + cddb);
            xcssb = (capbs + param.B2GSoverlapCap - (cgsb + cbsb + cdsb));
            xcggb = (cggb + param.B2GDoverlapCap
          + param.B2GSoverlapCap
          + param.B2GBoverlapCap);
            xcgdb = (cgdb - param.B2GDoverlapCap);
            xcgsb = (cgsb - param.B2GSoverlapCap);
            xcbgb = (cbgb - param.B2GBoverlapCap);
            xcbdb = (cbdb - capbd);
            xcbsb = (cbsb - capbs);

            state.Matrix[B2gNode, B2gNode] += xcggb * state.Omega;
            state.Matrix[B2bNode, B2bNode] += (-xcbgb - xcbdb - xcbsb) * state.Omega;
            state.Matrix[B2dNodePrime, B2dNodePrime] += xcddb * state.Omega;
            state.Matrix[B2sNodePrime, B2sNodePrime] += xcssb * state.Omega;
            state.Matrix[B2gNode, B2bNode] += (-xcggb - xcgdb - xcgsb) * state.Omega;
            state.Matrix[B2gNode, B2dNodePrime] += xcgdb * state.Omega;
            state.Matrix[B2gNode, B2sNodePrime] += xcgsb * state.Omega;
            state.Matrix[B2bNode, B2gNode] += xcbgb * state.Omega;
            state.Matrix[B2bNode, B2dNodePrime] += xcbdb * state.Omega;
            state.Matrix[B2bNode, B2sNodePrime] += xcbsb * state.Omega;
            state.Matrix[B2dNodePrime, B2gNode] += xcdgb * state.Omega;
            state.Matrix[B2dNodePrime, B2bNode] += (-xcdgb - xcddb - xcdsb) * state.Omega;
            state.Matrix[B2dNodePrime, B2sNodePrime] += xcdsb * state.Omega;
            state.Matrix[B2sNodePrime, B2gNode] += xcsgb * state.Omega;
            state.Matrix[B2sNodePrime, B2bNode] += (-xcsgb - xcsdb - xcssb) * state.Omega;
            state.Matrix[B2sNodePrime, B2dNodePrime] += xcsdb * state.Omega;
            state.Matrix[B2dNode, B2dNode] += gdpr;
            state.Matrix[B2sNode, B2sNode] += gspr;
            state.Matrix[B2bNode, B2bNode] += gbd + gbs;
            state.Matrix[B2dNodePrime, B2dNodePrime] += gdpr + gds + gbd + xrev * (gm + gmbs);
            state.Matrix[B2sNodePrime, B2sNodePrime] += gspr + gds + gbs + xnrm * (gm + gmbs);
            state.Matrix[B2dNode, B2dNodePrime] -= gdpr;
            state.Matrix[B2sNode, B2sNodePrime] -= gspr;
            state.Matrix[B2bNode, B2dNodePrime] -= gbd;
            state.Matrix[B2bNode, B2sNodePrime] -= gbs;
            state.Matrix[B2dNodePrime, B2dNode] -= gdpr;
            state.Matrix[B2dNodePrime, B2gNode] += (xnrm - xrev) * gm;
            state.Matrix[B2dNodePrime, B2bNode] += -gbd + (xnrm - xrev) * gmbs;
            state.Matrix[B2dNodePrime, B2sNodePrime] += -gds - xnrm * (gm + gmbs);
            state.Matrix[B2sNodePrime, B2gNode] += -(xnrm - xrev) * gm;
            state.Matrix[B2sNodePrime, B2sNode] -= gspr;
            state.Matrix[B2sNodePrime, B2bNode] += -gbs - (xnrm - xrev) * gmbs;
            state.Matrix[B2sNodePrime, B2dNodePrime] += -gds - xrev * (gm + gmbs);
        }
    }
}
